Discrete breathers in nonlinear magnetic metamaterials. 
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Magnetic metamaterials composed of split-ring resonators or U— type elements may exhibit dis- 
creteness effects in THz and optical frequencies due to weak coupling. We consider a model one- 
dimensional metamaterial formed by a discrete array of nonlinear split-ring resonators with each 
ring interacting with its nearest neighbours. On-site nonlinearity and weak coupling among the 
individual array elements result in the appearence of discrete breather excitations or intrinsic lo- 
calized modes, both in the energy-conserved and the dissipative system. We analyze discrete single 
and multibreather excitations, as well as a special breather configuration forming a magnetization 
domain wall and investigate their mobility and the magnetic properties their presence induces in 
the system. 

PACS numbers: 63.20.Pw, 75.30.Kz, 78.20.Ci 



C3 



-a 
o 

(N 
> 

in 
o 
\o 
o 

-I— > 
c3 



1 

C 

O 

o 



X 



Artificial non-magnetic materials exhibiting magnetic 
properties in the Terahertz and optical frequencies have 
been recently predicted theoretically^, 0] and demon- 
strated experimentally^ 0, 0. The key element for 
most of these magnetic metamaterials (MMs) is ei- 
ther the split-ring resonator (SRR) or its [/—shaped 
modification 6] . The realization of MMs at such (and 
possibly higher) frequencies will affect substantially THz 
optics and their applications in devices of compact cavi- 
ties, adaptive lenses, tunable mirrors, isolators and con- 
verters. Moreover, MMs with negative magnetic re- 
sponse can be combined with plasmonic wires that ex- 
hibit negative permittivity, producing left-handed mate- 
rials (LHM), i.e. metamaterials with negative magnetic 
permeability \i and dielectric permitiyity e leading to a 
negative index of refraction |7ll8ll3. llfHll| . In the present 
work we focus entirely on MMs that have the additional 
features of being nonlinear as well as discrete. While non- 
linearity results in self-focusing, discreteness induces lo- 
calization and, as a result, the combination of both leads 
in the generation of nonlinearly localized modes of the 
type of Discrete Breather (DB) (HE El 13 E3 These 
modes act like stable impurity modes that are dynami- 
cally generated and may alter propagation and emission 
properties of the system. 

We consider a planar one-dimensional (ID) array of 
N identical SRRs with their axes perpendicular to the 
plane; each unit is equivalent to an RLC oscillator with 
sclf-inductance L, Ohmic resistance R and capacitance 
C. The units become nonlinear due to the Kerr di- 
electric that fills their gap and has permittivity equal 
to e(|E| 2 ) = e (e e + a\'E\ 2 /E 2 ), where E is the elec- 
tric component of the applied electromagnetic field, E c 
is a characteristic electric field, q the linear permittiv- 
ity, eo the permittivity of the vacuum, and a = +1 



(a = —1) corresponding to self- focusing (-defocusing) 
nonlinearity 0, 0, 0. As a result, the SRRs aquire 
a field-dependent capacitance C(|E| 2 ) = e(|E g | 2 )^4/d 5 , 
where A is the area of the cross-section of the SRR wire, 
E g is the electric field induced along the SRR slit, and 
d g is the size of the slit. Since C(U n ) = dQ n /dU n , the 
charge Q n stored in the capacitor of the n-th SRR is 
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U n , n=l,2,...,N, 



(1) 



where Ci — eaepA/dg is the linear capacitance, U n — 
dgEg n is the voltage across the slit of the nth SRR, and 
U c = d g E c . Neighbouring SRRs are coupled due to mag- 
netic dipole-dipole interaction through their mutual in- 
ductance M, which decays as the cube of the distance. 
For weak coupling between SRRs in a planar configura- 
tion, it is a good approximation to consider only nearest 
neighboring SRR interactions. Then, the dynamics of Q n 
and the current /„ circulating in the nth SRR is described 
by 
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where £ is the electromotive force (emf ) induced in each 
SRR due to the applied field, and f(Q n ) — U n is given 
implicitly from Eq. Q). The value of £ at a given instant 
is proportional to the magnetic field component of the 
applied field perpendicular to the SRR plane, and/or the 
electric field component parallel to the side of the SRRs 
which contains the slit[2Cj|. Using the relations ujJ 2 = 
LCi, t = tup, I c = U c Lu e Ce, Q c = C e U c , £ = U c e, I n = 
I c ini Qn = Qc(In, Eqs. and 10 can be normalized to 

dq n _ . 



(4) 
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— (Xin-i - in + Ai n+ i) = ji n - e(r) + f{q n ), (5) 

where 7 = RCeOJe is the loss coefficient, and A = M/L is 
the coupling parameter. In the following, we use periodic 
boundary conditions (i.e., ijv+i = i\, *o = *jv) except 
otherwise stated. Analytical solution of Eq. for it„ = 
/(<?„) with the conditions of u n being real and u n (q n = 
0) = 0, gives the approximate expression 

f(q n ) ~ g„ - -% 3 + 3 f ' q 5 n + 0{q 7 n ) (6) 
dee \^ e f/ 

That is, the on-site potential V(q n ) = Jq" f{q' n )dq' n is 
soft for focusing nonlinearity and hard for defocusing 
nonlinearity. Substituting q n — A cos (kDn — f2r) into 
the linearized Eqs. and (JSJ with e = and 7 = 0, we 
obtain the frequency spectrum of linear excitations 

n k = [1-2A cos(fcL>)p 1/2 , (7) 

where £1 = lu/loi is the normalized frequency, D is the 
separation of neighbouring SRR centers (unit cell size), 
and k the wavenumber (— ir < k D < n). 

The parameter A can be calculated numerically for any 
SRR geometry, since the magnetic field of the current 
circulating the SRR is well known. Here we estimate 
A with a simple model 0, neglecting the effects of non- 
linearity and coupling on the resonance frequency |l8| . 
For not very small array dimensions, the inductance of 
a circular SRR of radius a with circular cross-section of 
diameter h, is L = fiQa[hi(l6a / h) — 1.5], where /i is 
the permeability of the vacuum. For a squared SRR 
with square cross-section with side length £ = 5 ixm, 
t = w = d g = 1 fim the SRR depth, width, and slit 
size, respectively, length of unit cell D = 7 ^m0, and 
using that £' = 4(£ — w) — d g is the length of the axis 
of the wire, a = £'/2ir, h = a/4 wt/w, we arrive at 
L ~ 6 x 10~ 12 H. For this L, we evaluate the capac- 
itance necessary for providing the resonance frequency 
for a single SRR, f r ~ l/2irs/LC~i = 6.2 THz, consis- 
tent with the available experimental informational, to 
be Ce — 11 x 10~ 17 F. Consider two neighbouring SRRs 
(1 and 2) in an array of circular SRRs of radius a with 
circular cross-section of diameter h. The flux $2 thread- 
ing SRR 2 due to the induced magnetic field in SRR 1 
Bi(r) ~ /x o <S7i/47rr 3 -|-0((a/r) 3 ), where Ii is the induced 
current in SRR 1, S = ira 2 is the SRR area and r is the 
distance from its center (r <~ D), is approximately $2 — 
Bi(r = D)S. Then, M = <f> 2 /h - HoS 2 /4ttD 3 and A ~ 
(n/4:)(a/D) 3 /[ln(l6a/h) — 1.5]. For an array of squared 
SRRs with square cross-section with dimensions as in 
we obtain A ~ (£' /L>) 3 /32vr 2 [?n(4f /Vmvt) - 1.5] ~ 0.02. 
For silver made SRRs, whose conductivity and skin depth 
are a ~ 6.15 x 10 7 S/m and 8 ~ 20 nm, respectively, we 
obtain i? = 2a/cr/i(5 = £' /2a5^/^i ~ 3.44, and 7 ~ 0.01. 




50 14 

t 



FIG. 1: (color online). Time evolution of a Hamiltonian 
breather for approximately two periods for A = 0.02, Tj, = 
6.69. a = +1. e/ = 2 and AT = 50. 




FIG. 2: Moving Hamiltonian breathers. Right: breather 
amplitude for T5 = 6.69, a = +1, ee = 2, N = 50, and 
A = 0.062 Left: space-time evolution of the center of energy 
Xe for the breather shown in the right figure, for A p = 0.1 
(solid line); 0.2 (dotted line); 0.3 (dashed line). 

We consider first the lossless case without applied field 
(7 = 0, e = 0). Then, Eqs. 10} and JSJ can be derived 
from the Hamiltonian 

n = Y,{\^n + V iln)-\qnqn+l\. (8) 

For Hamiltonian systems DBs may be constructed from 
the anti-continuous (AC) limit [14|, where all oscillators 
are uncoupled (A — * 0), obeying identical dynamical 
equations. Fixing the amplitude of one of them (say the 
one located at n = n^) to a specific value qi,, with the 
corresponding current % = 0, we can determine the oscil- 
lation period T^. An initial condition with q n = for any 
n ^ rib, q-n b = qbj an d q n — i n — for any n, represents a 
trivial DB. Continuation of this solution for A ^ using 
the Newton's method ^3], results in numerically exact 
DBs up to Xmax, where the linear excitation frequency 
band (which expands with increasing A) reaches the DB 
frequency u>b = 27r/T&. The linear stability of Hamil- 
tonian DBs is addressed through the eigenvalues of the 
monodromy matrix (Floquet coefficients). Fig. 1 shows 
the time evolution of a typical, linearly stable, highly lo- 
calized, DB excitation (X m ax ~ 0.067 for the chosen pa- 
rameters). In this figure, plotted vs. time t and array site 
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n, is the normalized current i n circulating the nth SRR. 
Another trivial DB can be obtained for q n = for any 
n 7^ nb, q nb = and q n = i n = for any n, correspond- 
ing to what we could call a "dark" DB in analogy with 
the dark soliton in nonlinear continuous systems. Such 
a DB can be continued up to A ~ 0.025 but it is linearly 
unstable except for very small A. In order to investigate 
the mobility of these DBs we followed the procedure de- 
scribed be Chen et al\lj^. According to this work, in 
order to generate a (steady state) moving DB, having 
obtained a static DB (q°,i° = 0) by Newton's method, 
we integrate Eqs. (@J and J5J using as initial condition 
(q(r = 0), i(r = 0)) = (q°, i° = 0) + A p (0, Si), where X p 
is the perturbation strength, and the perturbation vector 
Si corresponds to the current part of the (normalized to 
unity) pinning mode eigenvector. The resulting DB mo- 
tion is followed by plotting the instantaneous center of 
localization of energy Xe of the DB for several values of 
X p and a A value close to X m ax (Fig- 2); the parameter 
Xe is defined as 

N 

X E = J2 n - E n/ E tot, (9) 
n=l 

where E n is the energy at site n and E to t = Y] „_i E n . 
We note that Hamiltonian DBs move slowly through the 
lattice as a result of the perturbation. Their velocity de- 
creases with increasing A, although not uniformly with 
X p ; in particular, the DB velocity as a function of A de- 
creases faster as A p increases. 

In order to generate DBs for the forced and damped 
system we start by solving Eqs. |@J and JSJ in the AC 
limit with emf e(r) = £q sin(Jlr) . We identify two 
different amplitude attractors of the single SRR oscilla- 
tor, with amplitudes q n ~ 1.6086 and qi ~ 0.2866 for 
the high and low amplitude attractor, respectively. Sub- 
sequently, we fix the amplitude of one of the oscillators 
(say the one at n = rib) to q n and all the others to qi 
{i n are all set to zero). Using this configuration as ini- 
tial condition, we turn on adiabatically the coupling A. 

The initial condition can be continued for A ^ lead- 
ing to dissipative DB formation . The time evolution 
of a typical dissipative DB is shown in Fig. 3. Both 
the DB and the background are oscillating with different 
amplitudes (high and low, respectively). This should be 
compared to the Hamiltonian DB in Fig. 1, where the 
background is always zero. By interchanging q^ and qi in 
the initial conditions, we obtain another DB oscillating 
with low amplitude, while the background oscillates with 
high amplitude (Fig. 4). With appropriate initial con- 
ditions we can also obtain multi-breathers where two or 
more sites oscillate with high (low) amplitude, while the 
other ones with low (high) amplitude. Next, we fix the 
amplitude of half of the SRRs in the array (say those for 
n > N/2) to qh and the others to qe, and integrate Eqs. 
PJ) and JSJ from the AC limit with open-ended boundary 




FIG. 3: (color online). Time evolution of a one-site dissipa- 
tive breather during approximatelly two periods for Tj, = 6.82, 
A = 0.02, 7 = 0.01, £ = 0.04, a = +1, e t = 2 and N = 50. 
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FIG. 4: (color online). Time evolution of a one-site dissi- 
pative breather of the second type (see text) during approx- 
imatelly two periods, for A = 0.01 and the other parameters 
as in Fig. 3. 



conditions (i.e., ijv+i = io = 0). In this way we obtain 
an oscillating domain- wall, as shown in Fig. 5. 

The Hamiltonian DBs investigated are linearly stable, 
ensuring that they are not affected by small amplitude 
perturbations. On the other hand, the dissipative DBs 
are attractors for initial conditions in the correspond- 
ing basin of attraction and are robust against different 
kinds of small perturbations 15J. We analysed numeri- 
cally and confirmed the stability of the DBs presented 
above under various kinds of perturbations; we followed 
the purturbed DB evolution for long time intervals (over 
2 x 10 4 Th) without observing any significant change in 
the DB shape H3. 

The dissipative system, which includes forcing due to 
the applied field, offers the possibility to study its mag- 
netic response. Assume that the emf is induced by the 
magnetic component H = Hq cos(wt), of the applied 
field. Then, at least for uniform solutions (J„ = /), 
the magnetization M = SI/D 3 can be defined. In the 
direction perpendicular to the SRRs plane, the general 
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FIG. 5: (color online). Time evolution of a domain- wall 
breather during approximately two periods. Parameters as 
in Fig. 3. 




FIG. 6: (color online). Time evolution of Ki n {r) (red curve), 
compared with cos(f2r) (black curve), and their sum (green 
curve), for two SRRs of the domain- wall DB excitation, (a) 
low amplitude current oscillation SRR (n = 15); (b) high 
amplitude current oscillation SRR (n = 35). Parameters as 
in Fig. 3. 



and much larger in magnitude than that. Thus, for large 
enough SRR arrays, one may obtain domain-wall DBs 
connecting domains of the array with positive and nega- 
tive /i. 

In conclusion, a ID planar array of nonlinear SRRs 
coupled through nearest-neighbor mutual inductancies 
was investigated numerically. The existence of DBs of 
various types, for both the energy-conserved and the 
dissipative system, was demonstrated. We found that 
longer range interaction does not affect the DB properties 
substantially ji^. We found that Hamiltonian DBs may 
be set into uniform motion under a small perturbation. 
We also obtained a special DB solution (magnetization 
domain- wall) , which separates domains of the array with 
different magnetization. Multiple magnetization states 
are possible in this system due to nonlinearity, which al- 
lows either for negative or positive /i below resonance. 
Moreover, one can exploit multibreathers and domain- 
wall DBs to create MMs with domains of opposite sign 
magnetic response. Discreteness effects may appear in 
SRR arrays with dimensions close to those reported in 
, even though the field wavelength is much larger than 
the array dimensions. 

We acknowledge support from the grant "PYTHAGO- 
RAS II" (KA. 2102/TDY 25) of the Greek Ministry of 
Education and the European Union. 



relation B = /i (H + M) gives 

B = S (cos(r2r) + «i(r)), (10) 

where B = SqU c /SCIlj^, and k = fioS 2 ft/eoD 3 L. For the 
material parameters used above k ~ 3. From Eq. i|l(J|) 
negative magnetic response appears whenever the second 
term in the parentheses is larger in magnitude than the 
first one, and has opposite sign. Then, one may assign 
a negative fi to the medium. Without nonlinearity, the 
unique, uniform solution for the SRR array gives positive 
response below the resonance frequency (~ u>i). Nonlin- 
earity allows the existence of multiple stable states, which 
makes it possible to obtain either positive or negative \x 
below u>£, depending on the state of the system. More- 
over, exploiting DB excitations, MMs with domains of 
opposite sign magnetic responses can be created. In Fig. 
6 we show the time evolution of cos(fir) and «i n (r) as 
well as their sum, for two SRRs of the domain-wall DB 
(Fig. 5), relatively far from the domain- wall and the 
ends. The SRR with low amplitude current oscillation 
(ri = 15) shows positive (paramagnetic) response. In 
contrast, the SRR with high amplitude current oscilla- 
tion (n — 35) shows extreme diamagnetic (negative) re- 
sponse, since kj„(t) is almost out of phase with cos(fir) 
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